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By combining classical Monte Carlo and Bethe ansatz techniques we devise a numerical method to construct 
the Truncated Generalized Gibhs Ensemble (TGGE) for the spin-| isotropic Heisenberg (XXX) chain. The key 
idea is to sample the Hilhert space of the model with the appropriate GGE probability measure. The method can 
be extended to other integrable systems, such as the Lieh-Liniger model. We benchmark the approach focusing 
on GGE expectation values of several local observables. As finite-size effects decay exponentially with system 
size, moderately large chains are sufficient to extract thermodynamic quantities. The Monte Carlo results are 
in agreement with both the Thermodynamic Bethe Ansatz (TBA) and the Quantum Transfer Matrix approach 
(QTM). Remarkably, it is possible to extract in a simple way the steady-state Bethe-Gaudin-Takahashi (BGT) 
roots distributions, which encode complete information about the GGE expectation values in the thermodynamic 
limit. Einally, it is straightforward to simulate extensions of the GGE, in which, besides the local integral of 
motion (local charges), one includes arbitrary functions of the BGT roots. As an example, we include in the 
GGE the first non-trivial quasi-local integral of motion. 


Introduction .— The issue of how statistical ensembles 
arise from the out-of-equilibrium dynamics in isolated quan¬ 
tum many-body system is still a fundamental, yet challeng¬ 
ing, problem. The main motivation for the renewed inter¬ 
est in this topic is the high degree of control reached in 
out-of-equilibrium experiments with cold atomic gases*^'^. 
The paradigm experiment is the so-called global quantum 
quench^^, in which a system is initially prepared in an eigen¬ 
state Itho) of a many-body Hamiltonian %. Then a global 
parameter of Tf is suddenly changed, and the system evolves 
unitarily under the new Hamiltonian Tf'. At long times af¬ 
ter the quench the system reaches a steady state, as it has 
been confirmed by experiments^. In integrable models the 
presence of non-trivial local conserved quantities, besides the 
energy, strongly affects the dynamics and the nature of the 
steady state. As for now, despite the tremendous theoretical 
effort'^^^®, it is still unclear whether such steady-state can be 
described by a statistical ensemble, and how to construct it. 

It has been proposed that the long-time stationary value of a 
generic local operator O is described by a Generalized Gibbs 
Ensemble'^’^^ (GGE) as {O) = Here ex¬ 

tends the Gibbs density matrix by including all the extra con¬ 
served quantities Xj (charges) as 

pGGE = ^-1 exp ( - \jXj). (1) 

In (I), and in the rest of the paper, repeated indices are 
summed over. Z is a normalization factor. The Ay are La¬ 
grange multipliers to be fixed by imposing (T'o|Xy I'I'o) = 
{Xj), and X2 = %' is the post-quench Hamiltonian. In re¬ 
alistic situations one deals with the truncated GGE"'® (TGGE), 
i.e., considering only the “most local” charges. 

While the validity of the GGE has been largely confirmed 
in non-interacting theories-'®’^^’"'®’®^’®®, in interacting ones the 
scenario is far less clear (see Ref. 48 for numerical results in 
an interacting spin chain). Eor Bethe ansatz solvable mod¬ 
els the so-called Quench Action method"'"' allows for an ex¬ 
act treatment of the steady state, provided that the overlap 
between the initial state |T'o) and the eigenstates of Tf' are 
known. In several cases the Quench Action is in disagreement 


with the TGGE®'’®^“®®, whereas it is supported by numerical 
simulations®^. The origin of this discrepancy remained un¬ 
known until very recently. In Ref. 56 it has been shown that it 
is possible to “repair” the GGE by including the quasi-Iocal 
charges®®”®*. Remarkably, this repaired GGE is in perfect 
agreement with the Quench Action®®, confirming that the de¬ 
scription of the steady state with the GGE is correct, provided 
that the appropriate set of local and quasi-local charges is con¬ 
sidered. 

On the other hand, numerical methods, such as the time de¬ 
pendent density matrix renormalization group®^ ®® (tDMRG), 
have been mostly used to simulate the post-quench dynam¬ 
ics in microscopic models. However, no numerical attempt 
to explore the GGE per se has been undertaken yet. The 
aim of this work is to provide a Monte-Carlo-based frame¬ 
work for studying the GGE, and its possible extensions, in 
Bethe ansatz solvable models. We restrict ourselves to finite- 
size systems. Thermodynamic quantities can be extracted by 
a standard finite-size scaling analysis. Moreover, as finite- 
size corrections decay exponentially with system size®"', mod¬ 
erately large systems are sufficient to access the thermody¬ 
namic limit. The method relies on the detailed knowledge 
of the Hilbert space structure provided by the Bethe ansatz 
formalism, and on the Bethe-Gaudin-Takahashi (BGT) equa¬ 
tions®®’®®. The key idea is to sample the model Hilbert space 
according to the GGE probability measure given in (1). We 
should mention that the same idea has been already explored 
in Ref. 67 for the Gibbs ensemble. The method allows one 
to obtain GGE expectation value for generic observables, pro¬ 
vided that their expression in terms of the roots of the BGT 
equations are known. Remarkably, it is also possible to ex¬ 
tract the steady-state roots distributions, which encode the 
complete information about the (GGE) ensemble. It is also 
straightforward to extend the GGE including in (1) arbitrary 
functions of the BGT roots. This could be useful, for in¬ 
stance, to investigate the effects of quasi-local charges. Ei¬ 
nally, we should mention that, in principle, GGE averages of 
local observables can be computed using exact diagonaliza- 
tion or Quantum Monte Carlo. However, both these methods 
require the operatorial expression of the conserved charges 
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FIG. 1. The Generalized Gibbs Ensemble (GGE) for the Heisenberg spin chain with L = 16 sites: Numerical results obtained using the 
Hilbert space Monte Carlo sampling approach. Only the first two even conserved charges 12 , Tr and the first quasi-local one Tqs are included 
in the GGE. T 2 is the Hamiltonian. In all the panels different symbols correspond to different values of the Lagrange multipliers A 4 , Xqs- The 
circles correspond to the Gibbs ensemble, i.e., A 4 = Xqs = 0. The a:-axis shows the inverse temperature A 2 = (a) The GGE average 

{X 2 IL). (b) Variance of the GGE fluctuations a^{X 2 )/L = ((TI) ~ (T 2 )^)/T as a function of /3. (c)(d) and (e)(f): Same as in (a)(b) for 
Xi and Te, respectively. In all panels the lines are the Quantum Transfer Matrix (QTM) results, (g) x!P plotted versus /3, with x being the 
magnetic susceptibility per site. 


(see Ref. 68 for the XXX chain), whereas our results rely 
only on their expression (typically simple) in terms of the 
BGT roots. 

We benchmark the approach focusing on the spin-1 
isotropic Heisenberg chain {XXX chain), which is the ven¬ 
erable prototype of integrable models®®. We consider several 
TGGEs (cf. (1)) constructed including X 2 ,l 4 , and the first of 
the recently discovered®®’®' quasi-local charges Xqs {Hi in 
Ref. 56). We focus on the conserved charges averages {Xj/L), 
and on their ensemble fluctuations a^{Xj) = (2|) — 
which are related to well-known physical observables, such 
as the energy density, and the specific heat. We also compute 
the spin susceptibility per site x- Already for a chain with 
L — sites the Monte Carlo data perfectly agree with both 
the standard Thermodynamic Bethe Ansatz®® (TBA) and the 
Quantum Transfer Matrix approach"'®’®' (QTM). Notice that 
this is the first direct numerical verification of the QTM ap¬ 
proach in the XXX chain. Finally, we extract the BGT roots 
distributions for both the Gibbs ensemble and the GGE. In 
both cases the finite-size effects are negligible for small roots, 
which are the relevant ones to describe the long-wavelength 
physics. For the Gibbs ensemble we compare our numerical 
data with standard finite-temperature Thermodynamic Bethe 
Ansatz (TBA) results, finding excellent agreement. 

The-Heisenberg-spin-chain .— The XXX chain with L 
sites is defined by the Hamiltonian 




i=l 


^{S+S-^, + S-St,,) + StS!+, 



where Sf = (erf ± iaf)/2 are spin operators acting on the 
site i, Sf = erf /2, and erf the Pauli matrices. We fix J = 1 
and use periodic boundary conditions, identifying sites L + 1 
and 1. The total magnetization S'f = — M, 

with M number of down spins (particles), commutes with (2), 
and it is here used to label its eigenstates. 

In the Bethe ansatz formalism each eigenstate of (2) is uni- 
vocally identified by M parameters {xa G 'C}a=i- '^he 


limit L —> c» they form “string” patterns along the imaginary 
axis of the complex plane (string hypothesis®®’®®). Strings of 
length 1 < n < M (so-called n-strings) are parametrized as 
= Xn-'j — i{n — 1 — 2j). Here G M is the string 
real part (string center), j = 0,1,..., n — 1 labels different 
string components, and 7 denotes different string centers. The 
string hypothesis is not correct for finite chains, although de¬ 
viations typically decay exponentially with L. Physically, the 
n-strings correspond to eigenstate components containing n- 
particle bound states. The {xn-.y} are obtained as the roots of 
the Bethe-Gaudin-Takahashi (BGT) equations®® ®® 

XP Yi{x — 27r/72.-y-p ^ ^ Qm,n (■t^n;7 (3) 

(m,/3)^(n,j) 


Here = 2arctan(a:/n), Qm,n{x) is the scattering 

phase between different roots®®, and In-^ G are the so- 
called Bethe-Gaudin-Takahashi quantum numbers. The 
satisfy the upper bound \In-,-y \ < ImAxin, L, M), with /max a 
known function®® of n, L, M. Every choice of /„;-y identifies 
an eigenstate of (2). We define the “string content” of each 
eigenstate as 5 = {si,..., sm}, with 0 < s„ < [M/nJ the 
number of n-strings. The local conserved charges Xj of the 
XXX chain are given as 


2(j+i = 




(j - 1)! dy3 


logA(y) 


y^i 


( 4 ) 


where A(j/) is the eigenvalue of the quantum transfer matrix®®, 
with y a spectral parameter. X2 is the XXX Hamiltonian. 
The analytic expression of Xj in terms of the Pauli matrices 
is known®* for j < 10. The support of Xj, i.e., the num¬ 
ber of adjacent sites where Xj acts non trivially, increases lin¬ 
early with j, i.e., larger j correspond to less local charges. 
The eigenvalues of Xj on a generic eigenstate are obtained by 
summing the contributions of the different BGT roots inde¬ 
pendently. For instance, the energy eigenvalue is obtained as 
E — 2 7 xiliv? -F A similar result holds true for 

the quasi-local charges®®. 
















3 


The-Hilbert-space-Monte-Carlo-sampling .— For a finite 
chain the GGE (cf. (1)) can be obtained by importance sam¬ 
plings^ of the eigenstates of (2). One starts with an initial 
M-particle eigenstate, with string content S = {si, ..., sm }, 
and identified by a BGT quantum number configuration C = 
{In-,-y}n=i (7 = 1, ■ ■ •, s„). The Corresponding charges 
eigenvalues are {2y}. Then a new eigenstate is generated with 
a Monte Carlo scheme. Each Monte Carlo step (mcs) consists 
of three moves: 


Choose a new particle number sector M' , and string 
content S' = {s^,..., s'^/} with probability S'* 

r{M',s') 


r{M',s') = 


B{L,LI2) 


M' 

X{b{U 


s') 


(5) 


i=l 


Generate a new quantum number configuration C com¬ 
patible with the S' obtained in step 1. Solve the corre¬ 
sponding BGT equations (3). 


Calculate the charge eigenvalues X' and accept the new 
eigenstate with the Metropolis probability: 



1 , 


L - 2M' + 1 
L-2M + I 




}■ 


( 6 ) 


In (5) B{x, y) = x\/{y\{x — t/)!) is the Newton binomial and 

Ci = L — J2jLi tijB'j, with Uj = 2Min(i, j) — S^j. In ( 6 ) the 
factor in front of the exponential takes into account that Xj and 
the observables that we consider are invariant under SU{2) ro¬ 
tations. Crucially, the steps 1 and 2 are necessary to account 
for the density of states of the model (equivalently, the Yang- 
Yang entropy, see below), and are the same as for the Gibbs 
ensemble^s xhe iteration of 1-3 defines a Markov chain, 
which, after some thermalization steps, generates eigenstates 
distributed according to (1). Interestingly, by trivially modi¬ 
fying ( 6 ) it is possible to simulate more exotic ensembles in 
which, in addition to Xj, one considers arbitrary functions of 
the BGT roots. The GGE average {O) of a generic operator is 
obtained as 


{O) 


lim 


1 


tvmcs ^ OO JYm 


|s> 


(7) 


where iVmcs is the total number of eigenstates |s) sampled in 
the Monte Carlo. Moreover, for all the observables consid¬ 
ered here the contributions of the BGT roots can be summed 
independently, i.e., 

{s\0\s) =^foix„-j), ( 8 ) 

n,7 
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FIG. 2. Finite-size scaling of the GGE averages in the Heisenberg 
chain: Numerical results obtained from the Hilbert space Monte 
Carlo sampling. Here the GGE is constructed including X 2 ,77, with 
Lagrange multipliers A 2 = /3, A 4 = 1. (a) {X 2 / 7 /) plotted versus 
the chain size L for several values of /3. The dash-dotted lines are 
exponential fits, (b) Same as in (a) for 77. 


{{Xj) - (lj)2)/E (panels (b)(d)(f)). Panel (g) plots x//3, with 
X the spin susceptibility. Notice that {X 2 IL) is the energy 
density, while (t^(Z 2 )/E is related to the specific heat. In all 
panels the data correspond to the TGGE constructed with the 
first two even charges X2,X4, and the first non-trivial quasi¬ 
local charge Different symbols correspond to dif¬ 

ferent values of the associated Lagrange multipliers, namely 
A 4 = Xqs = 0 (Gibbs ensemble, circles in the Eigure), 
A 4 = 1 and Xqs = 0 (squares), and A 4 = 0 , Xqs = 1 
(crosses). In all panels the x-axis shows the inverse tem¬ 
perature X 2 = P- The data are Monte Carlo averages with 
Nmcs = 5 • 10® (cf. (7)). As expected, the different ensembles 
give different expectation values, implying that local observ¬ 
ables are able to distinguish different GGEs. In Eig. 1 the 
continuous lines are the analytic results obtained in the ther¬ 
modynamic limit using the QTM approach. These fully match 
the Monte Carlo data, signaling that finite-size effects are neg¬ 
ligible already for L = 16. 

The finite-size corrections are more carefully investigated 
in Eig. 2, plotting (I 2 ) and {X 4 ) (panels (a) and (b), respec¬ 
tively) versus /3. We focus on the TGGE with X2 = Xqs = 0 
and A 4 = 1. Clearly, finite-size effects decay exponen¬ 
tially®* with L for any /3. In (a) the dashed lines are fits to 
Cl -I- C 2 exp(—C 3 L), with Cl, C 2 , C 3 fitting parameters. Einite- 
size corrections are larger at lower temperature, and increase 
with the range of the operator (compare panels (a) and (b) in 
Eig. 2), as expected. 

Extracting-the-root-distributions. — In the thermody¬ 
namic limit in each n-string sector the roots of (3) become 
dense. Thus, instead of the eigenstates, one considers the 
corresponding root distributions p = {pn}'^^i- Eormally, 
Pn = \imL->.oo[L{xn-^+i “ Xn-j)]^^- The GGE average of a 
generic observable O becomes a functional integral as^°75 


where Xn-^ are the roots identifying the eigenstate |s), and 
fo{x) depends on the observable. 

The-GGE-for-local-observables .— The correctness of 
the Monte Carlo approach is illustrated in Eig. 1, considering 
the charge densities {XjjL) (panels (a)(c)(e) in the Eigure), 
and the variance of their ensemble fluctuations {Xj)/L = 


Tr{ exp (AjX,)C>} J Vpexp {^S[p] +XjXj\p])0[p]. (9) 

Here S[p] is the Yang-Yang entropy, which counts the number 
of eigenstates leading to the same p in the thermodynamic 
limit, and it is extensive. In (9) it is assumed that O becomes 
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FIG. 3. The root distributions p„ (x) (for n = 1,2,3) for the infinite temperature Gibbs (panels (a)-(c)) and the GGE equilibrium states (panels 
(d)-(f)): Numerical results for the Heisenberg spin chain obtained using the Hilbert space Monte Carlo sampling. Here the GGE is constructed 
including only X 2 and X 4 with fixed Lagrange multipliers A 2 = 0 and A 4 = 1. In all the panels the data are the histograms of the n-strings 
roots sampled in the Monte Carlo. The width of the histogram bins is Ax = 2/L, with L the chain size. In each panel different histograms 
correspond to different L. All the data are divided by 10^ for convenience. In (b) the arrow highlights the finite-size effects. In (a)-(c) the 
lines are the Thermodynamic Bethe Ansatz (TBA) results, (g) Einite-temperature effects: Monte Carlo data for for different values of 
the inverse temperature /?. 


a smooth functional of p in the thermodynamic limit. Eq. (8) 
becomes 

(5|0|s) ^ f dxpnix)fo{x). (10) 

n 

Since both S'!/)] and Xj[p] are extensive, the functional in¬ 
tegral in (9) is dominated by the saddle point^° p^^, with 
S{S + Xjlj)/Sp\p=psp = 0. Here p^P acts as a representa¬ 
tive state for the ensemble, and it contains the full information 
about the GGE equilibrium steady state. Eq. (7) and (10) im¬ 
ply that in the thermodynamic limit the histograms of the BGT 
roots sampled in the Monte Carlo converge to p^P. 

This is numerically supported in Eig. 3. Panels (a)-(c) plot 
the root distributions for n = 1, 2, 3 as a function of 

X for the representative state (saddle point) of the infinite- 
temperature Gibbs ensemble. In each panel the different his¬ 
tograms correspond to different chain sizes 18 < L < 30. 
The data are obtained using 5 • 10® Monte Carlo steps. The 
width of the histogram bins is varied with L as 2/L. In all 
the panels the full lines are the analytical Thermodynamic 
Bethe Ansatz®® (TBA) results. Clearly, deviations from the 
TBA vanish upon increasing the chain size (see for instance 
the arrow in panel (b)). Moreover, the corrections are larger 
on the tails of the distributions. This is expected since large 
roots correspond to large quasi-momenta, which are more sen¬ 
sitive to the lattice effects. Einally, finite-size effects increase 
with n, i.e., with the bound state sizes, as expected. The re¬ 
sults for the finite-temperature Gibbs ensemble are reported in 
Eig. 1 (g), for f3 = 1/2 and /3 = 1 (the different histograms). 
We focus on restricting ourselves to L = 30. The in¬ 

finite temperature histogram is reported for comparison. The 
continuous lines are now finite-temperature TBA results, and 
perfectly agree with the Monte Carlo data. Upon lowering 
the temperature the height of the peak at a; = 0 increases. 
This reflects that at /3 = 00 the tail of the root distributions 
vanish exponentially, whereas for /3 = 0 they are®® ^ 1/x^. 
Einally, panels (d)-(f) plot pn{x) for the TGGE constructed 


with 12,1-4 at fixed A 2 = 0, A 4 = 1 and for L = 30. In¬ 
terestingly, in contrast with the thermal case (see (a)), 
exhibits a double peak at small x. Similar to the infinite- 
temperature Gibbs ensemble ((a)-(c) in the Eigure), the data 
suggest that for L — 30 finite-size effects are negligible, at 
least for —2 < a; < 2 . 

Conclusions .— We presented a Monte-Carlo-based 
scheme for simulating the truncated Generalized Gibbs en¬ 
semble (TGGE) in finite-size integrable models. The key idea 
is the importance sampling of the model eigenstates using the 
GGE probability measure. The method relies on the Bethe 
ansatz formalism, and, in particular, on the Bethe-Gaudin- 
Takahashi (BGT) equations. The thermodynamic limit can 
be accessed by standard finite-size scaling analysis. Eor local 
quantities we observed that the finite-size corrections decay 
exponentially with the system size. Remarkably, the method 
allows to extract in a simple way the steady-state BGT root 
distributions, which contain full information about the (GGE) 
ensemble averages in the thermodynamic limit. Einally, it is 
possible to simulate extensions of the GGE, in which, besides 
the integral of motion, one includes arbitrary functions of the 
BGT roots. We benchmarked the method focusing on the 
spin-i isotropic Heisenberg chain. Specifically, we compared 
the Monte Carlo results with the standard Thermodynamic 
Bethe ansatz and the Quantum Transfer Matrix approach, 
finding excellent agreement. Einally, we simulated an 
extended GGE obtained by including the first non-trivial 
quasi-local charge. 

As an interesting research direction, we mention that it 
would be useful to generalize the method to simulate the GGE 
at fixed value of the conserved charges. This should be pos¬ 
sible using the standard microcanonical Monte Carlo tech¬ 
niques that have been developed in lattice gauge theory^® and 
in molecular dynamics simulations^^. Einally, by including 
in ( 6 ) the overlap contribution log |('I'o|^'j)|^ with Ifllo) the 
pre-quench initial state, and ) the eigenstates of the model, 
it should be possible to simulate the Quench Action^*^. 
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